# ------------------------------------- #
# Replication code for:
#
# Pomeroy, Caleb "Hawks Become Us: The Sense of Power and Militant Foreign Policy Attitudes," Security Studies.
#
# This script reproduces the supplementary analyses associated with the paper.
# ------------------------------------- #

# --- set your working directory --- #
setwd("~/Downloads/hawks_replication/")

# --- libraries --- #
library(readr)
library(ggplot2)
library(texreg)
library(mediation)
library(plyr)
library(scales)
library(caret)
library(ggeffects)
library(ggh4x)
library(mice)


# --- import data --- #
# for more details about variables included in the datasets, see "hawks_main_script.R" and the online appendix

# --- app data
app_survey <- read_csv("data/app_survey.csv")
app_survey$gender <- factor(app_survey$gender, levels = c("non_male", "male"))

# --- russia data
russia_survey <- read_csv("data/russia_survey.csv")
russia_survey$gender <- factor(russia_survey$gender, levels = c("female", "male"))
russia_survey$education <- factor(russia_survey$education, levels = c("bachelors", "less_bachelors"))

# --- qualtrics US data (survey #1, which included in-depth scenario)
qualtrics_survey_scenario <- read_csv("data/us_survey1.csv")
qualtrics_survey_scenario$gender <- factor(qualtrics_survey_scenario$gender, levels = c("non_male", "male"))

# --- china data
china_survey <- read_csv("data/china_survey.csv")
china_survey$gender <- factor(china_survey$gender, levels = c("female", "male"))

# --- qualtrics US data (survey #2, which included multiple measures of individual differences)
qualtrics_survey_inddiffs <- read_csv("data/us_survey2.csv")
qualtrics_survey_inddiffs$gender <- factor(qualtrics_survey_inddiffs$gender, c("non_male", "male"))
qualtrics_survey_inddiffs$education <- factor(qualtrics_survey_inddiffs$education, levels = c("bachelors", "less_bachelors"))


# --- Who feels that their state is powerful? --- #
# The big takeaway here is the lack of systematic patterns in predictors of a sense of state power,
# which is to say that no single subpopulation seems to drive the results.
summary(power_us <- lm(sense_power ~ nat_attach + education + gender + age +
                         harm_care + fairness_recip + ingroup_loyalty + authority_respect +
                         extraversion + agreeableness + conscientiousness + neuroticism + openness +
                         rwa + sdo, data = qualtrics_survey_inddiffs))

summary(power_china <- lm(sense_power ~ nat_attach +  education + gender + age +
                            harm_care + fairness_recip + ingroup_loyalty + authority_respect, data = china_survey))

summary(power_russia <- lm(sense_power ~ nat_attach + education + gender + age +
                             harm_care + fairness_recip + ingroup_loyalty + authority_respect, data = russia_survey))

# --- Table A2
screenreg(list(power_us, power_china, power_russia),
       booktabs = T,
       caption.above = T,
       caption = "OLS Results: Who Feels That Their State is Powerful?",
       custom.coef.names = c("(Intercept)",
                             "National Attachment",
                             "No College Degree",
                             "Male",
                             "Age",
                             "Harm/Care",
                             "Fairness/Reciprocity",
                             "Ingroup/Loyalty",
                             "Authority/Respect",
                             "Extraversion",
                             "Agreeableness",
                             "Conscientiousness",
                             "Neuroticism",
                             "Openness",
                             "RWA",
                             "SDO"),
       stars = c(.001, .01, .05, .10),
       symbol = "\\wedge",
       label = "table:ols_power_dv",
       custom.model.names = c("US Public",
                              "Chinese Public",
                              "Russian Public"))


# --- Factor Analysis: Sense of Power versus Militant Internationalism --- #
# Given that the paper expects the feeling of power to explain hawkishness, it is important to assess whether 
# state power is distinguishable from hawkishness (that is, whether these items load onto separate factors)
power_bind <- 
  cbind(qualtrics_survey_inddiffs$sense_power_great,
        qualtrics_survey_inddiffs$sense_power_wants,
        qualtrics_survey_inddiffs$sense_power_military,
        qualtrics_survey_inddiffs$sense_power_economic,
        qualtrics_survey_inddiffs$mi_demonstrate_resolve,
        qualtrics_survey_inddiffs$mi_strength_peace,
        qualtrics_survey_inddiffs$mi_war_unfortunate)
psych::fa.parallel(power_bind) # clearly two factors, one MI and one sense of power
factor.fit <- psych::fa(power_bind,nfactors=2, rotate="varimax",scores=TRUE, fm="pa")
# --- Table A3:
factor.fit$loadings 
# in short, the sense of power is the sense of capacity; hawkishness is what you do with that capacity


# --- Further Replications of Power’s Effect on MI --- #
# Figure 1 of the main text displays the relationship between the sense of power and MI for the three surveys
# that had various measures of individual differences (e.g., moral foundations, big 5).
# The models below simply show the same positive relationship between the sense of power and MI for the 
# APP survey and the Qualtrics survey that included the Iran scenario. 
summary(mi_app_simple <- lm(mi_war ~ sense_power, data = app_survey))
summary(mi_app <- lm(mi_war ~ sense_power + 
                       ideology + gender + age + urban_rural + education + white, data = app_survey))

summary(mi_qualtrics_simple <- lm(mil_int ~ sense_power, data = qualtrics_survey_scenario))
summary(mi_qualtrics_scenario <- lm(mil_int ~ sense_power + 
                                      ideology + gender + age + white + pid + nat_attach, data = qualtrics_survey_scenario))

# --- Table A5:
screenreg(list(mi_app_simple, mi_app, mi_qualtrics_simple, mi_qualtrics_scenario),
       booktabs = T,
       caption.above = T,
       caption = "OLS Results: Sense of Power Explains Militant Orientation (Further Replications)",
       custom.coef.names = c("(Intercept)",
                             "Sense of State Power",
                             "Conservative",
                             "Male",
                             "Age",
                             "Urban Resident",
                             "Less than College",
                             "White",
                             "Republican",
                             "National Attachment"),
       stars = c(.001, .01, .05, .10),
       symbol = "\\wedge",
       custom.model.names = c("(1)",
                              "(2)",
                              "(3)",
                              "(4)"))

# --- National Attachment Robustness Checks --- #
# For the correlational results, one concern is that perhaps individuals highly attached to their state are more likely 
# to positively evaluate their state as powerful and negatively evaluate others as threatening. Of course, the paper's models
# often control for national attachment. But, another way to think about this question is to simply subset individuals below 
# the median level of national attachment to check whether the same positive relationship exists between sense of power and MI.
summary(mi_us_simple <- lm(mil_int ~ sense_power, data = subset(qualtrics_survey_inddiffs, nat_attach < median(na.omit(qualtrics_survey_inddiffs$nat_attach)))))
summary(mi_us <- lm(mil_int ~ sense_power + 
                      nat_attach + education + gender + age + 
                      harm_care + fairness_recip + ingroup_loyalty + authority_respect +
                      extraversion + agreeableness + conscientiousness + neuroticism + openness +
                      rwa + sdo, data = subset(qualtrics_survey_inddiffs, nat_attach < median(na.omit(qualtrics_survey_inddiffs$nat_attach)))))

summary(mi_china_simple <- lm(mi_war ~ sense_power, data = subset(china_survey, nat_attach < median(na.omit(china_survey$nat_attach)))))
summary(mi_china <- lm(mi_war ~ sense_power + 
                         nat_attach +  education + gender + age + harm_care + fairness_recip + ingroup_loyalty + authority_respect, data = subset(china_survey, nat_attach < median(na.omit(china_survey$nat_attach)))))

summary(mi_russia_simple <- lm(mil_int ~ sense_power, data = subset(russia_survey, nat_attach < median(na.omit(russia_survey$nat_attach)))))
summary(mi_russia <- lm(mil_int ~ sense_power + 
                          nat_attach + education + gender + age + harm_care + fairness_recip + ingroup_loyalty + authority_respect, data = subset(russia_survey, nat_attach < median(na.omit(russia_survey$nat_attach)))))

# --- Table A6:
screenreg(list(mi_us_simple, mi_us, mi_china_simple, mi_china, mi_russia_simple, mi_russia),
       booktabs = T,
       caption.above = T,
       caption = "OLS Results: Sense of Power Explains Militant Orientation (Low National Attachment)",
       custom.coef.names = c("(Intercept)",
                             "Sense of State Power",
                             "National Attachment",
                             "Less than College",
                             "Male",
                             "Age",
                             "Harm/Care",
                             "Fairness/Reciprocity",
                             "Ingroup/Loyalty",
                             "Authority/Respect",
                             "Extraversion",
                             "Agreeableness",
                             "Conscientiousness",
                             "Neuroticism",
                             "Openness",
                             "Right-Wing Authoritarianism",
                             "Social Dominance Orientation"),
       stars = c(.001, .01, .05, .10),
       symbol = "\\wedge",
       custom.model.names = c("(1)",
                              "(2)",
                              "(3)",
                              "(4)",
                              "(5)",
                              "(6)"))


# --- Waltz Cross-Validated: Structure Still Explains A Lot With A Little --- #
# Akin to Waltz's structural intuitions, if we knew nothing about an individual except for their sense of state power, 
# how much variance in hawkishness could we predict? The cross-validation exercises below compare the fully-specified regressions
# against a simple regression using the sense of power as the only rhs predictor

# --- US public
set.seed(1912) 
train_control <- trainControl(method = "cv", number = 10)
ols_model_allvars_us <- train(mil_int ~., data =  na.omit(qualtrics_survey_inddiffs[,c("mil_int", "sense_power", "harm_care",
                                                                                    "fairness_recip", "ingroup_loyalty", "authority_respect",
                                                                                    "nat_attach", "rwa", "sdo",
                                                                                    "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness",
                                                                                    "education", "gender", "age")]), method = "lm",
                           trControl = train_control)

print(ols_model_allvars_us)

set.seed(1912) 
ols_model_onlypower_us <- train(mil_int ~., data =  na.omit(qualtrics_survey_inddiffs[,c("mil_int", "sense_power")]), method = "lm",
                             trControl = train_control)
print(ols_model_onlypower_us)

ols_model_onlypower_us$resample$model <- "only_power"
ols_model_onlypower_us$resample$survey <- "American Public"
ols_model_allvars_us$resample$model <- "all_vars"
ols_model_allvars_us$resample$survey <- "American Public"

# --- Chinese public
set.seed(1912) 
ols_model_allvars_china <- train(mi_war ~., data =  na.omit(china_survey[,c("mi_war", "sense_power", "harm_care",
                                                                            "fairness_recip", "ingroup_loyalty", "authority_respect",
                                                                            "nat_attach",
                                                                            "education", "gender", "age")]), method = "lm",
                                 trControl = train_control)
print(ols_model_allvars_china)

set.seed(1912) 
ols_model_onlypower_china <- train(mi_war ~., data =  na.omit(china_survey[,c("mi_war", "sense_power")]), method = "lm",
                                   trControl = train_control)
print(ols_model_onlypower_china)

ols_model_onlypower_china$resample$model <- "only_power"
ols_model_onlypower_china$resample$survey <- "Chinese Public"
ols_model_allvars_china$resample$model <- "all_vars"
ols_model_allvars_china$resample$survey <- "Chinese Public"

# --- Russian public
set.seed(1912) 
ols_model_allvars_russia <- train(mil_int ~., data =  na.omit(russia_survey[,c("mil_int", "sense_power", 
                                                                            "nat_attach",
                                                                            "harm_care", "fairness_recip", "ingroup_loyalty", "authority_respect",
                                                                            "gender", "age")]), method = "lm",
                                  trControl = train_control)
print(ols_model_allvars_russia)

set.seed(1912) 
ols_model_onlypower_russia <- train(mil_int ~., data =  na.omit(russia_survey[,c("mil_int", "sense_power")]), method = "lm",
                                    trControl = train_control)
print(ols_model_onlypower_russia)

ols_model_onlypower_russia$resample$model <- "only_power"
ols_model_onlypower_russia$resample$survey <- "Russian Public"
ols_model_allvars_russia$resample$model <- "all_vars"
ols_model_allvars_russia$resample$survey <- "Russian Public"

plot_df <- rbind(ols_model_onlypower_us$resample, ols_model_allvars_us$resample,
                 ols_model_onlypower_china$resample, ols_model_allvars_china$resample,
                 ols_model_onlypower_russia$resample, ols_model_allvars_russia$resample)

plot_df$model <- factor(plot_df$model, levels = c("all_vars", "only_power"))
levels(plot_df$model) <- c("All\nVariables", "Sense of\nPower Alone")
plot_df$survey <- factor(plot_df$survey, levels = c( "American Public", "Chinese Public", "Russian Public"))

# --- Figure A1
ggplot(plot_df, aes(x = model, y = RMSE, fill = model)) +
  geom_boxplot(color = "black", size = 1, width = .6, fatten = 1) +
  labs(y = "Root Mean Square Error\n(Smaller is Better)", x = NULL) +
  theme_bw() +
  facet_grid(~survey) +
  coord_cartesian(ylim = c(0,.3)) +
  scale_fill_manual(values = c("#99998e", "#635d67")) + 
  theme(axis.title.y = element_text(size = 18, margin = margin(r=15)),
        axis.text.y = element_text(size = 17),
        axis.text.x = element_text(size = 15),
        strip.text = element_text(size = 17, color = "gray94"),
        legend.position = "none",
        panel.spacing = unit(1, "lines"),
        panel.background = element_rect(fill = "#f6f6f5"),
        strip.background = element_rect(fill = "#54544c", color = "black", size = 1.5),
        plot.margin = margin(r=15, l=15, t= 15))

t.test(ols_model_onlypower_us$resample$RMSE, ols_model_allvars_us$resample$RMSE)
mean(ols_model_allvars_us$resample$RMSE)/mean(ols_model_onlypower_us$resample$RMSE)
t.test(ols_model_onlypower_china$resample$RMSE, ols_model_allvars_china$resample$RMSE)
mean(ols_model_allvars_china$resample$RMSE)/mean(ols_model_onlypower_china$resample$RMSE)
t.test(ols_model_onlypower_russia$resample$RMSE, ols_model_allvars_russia$resample$RMSE)
mean(ols_model_allvars_russia$resample$RMSE)/mean(ols_model_onlypower_russia$resample$RMSE)



# --- Re-Analysis of the 2019 Chicago Council–Levada Center Survey --- #
# Concerned about the small sample size of the Russia survey, this section re-analyzes a previous survey in Russia
# that included questions that can serve as a proxy for the sense of power, alongside various DVs that resemble hawkishness

# load data
chicago_levada_survey <- read_csv("data/chicago_levada_data.csv")
chicago_levada_survey$voted_putin <- factor(chicago_levada_survey$voted_putin)
chicago_levada_survey$gender <- factor(chicago_levada_survey$gender, levels = c("non_male", "male"))
# See the paper's online appendix as well as the original Chicago Council report for more info on the variables 

# create a sense of power proxy based on items that roughly capture this paper's concept of the sense of power
sense_russian_power <- chicago_levada_survey[,c("russia_rising", "russian_defensive_power", "russian_economy_improved", "russian_influence_improved")]
psych::fa.parallel(sense_russian_power) 
factor.fit <- psych::fa(sense_russian_power,nfactors=1, rotate="varimax",scores=TRUE, fm="pa")
factor.fit$loadings 
# basically the "Russia rising" item isn't a great fit with the other three items, but the results below are unchanged (if not slightly improved)
# when omitting that item. So, will keep all four items for theoretical reasons.
chicago_levada_survey$sense_power_factor <- factor.fit$scores[,1]
chicago_levada_survey$sense_power_factor <- scale(chicago_levada_survey$sense_power_factor)[,1]

# --- There's enough missing data to warrant multiple imputation
data_trim <- as.data.frame(chicago_levada_survey[,c("active_posture", "syria_involvement_good", "annexation_crimea_good",
                                            "luhansk_donetsk_part_russia", "sense_power_factor",
                                            "voted_putin", "gender", "age", "education", "income", "rural_region")],
                           stringsAsFactors = F)

imputed_data <- mice(data_trim, m=5, maxit = 50, method = 'pmm', seed = 1912)

# --- Now run the various regressions, pooling the estimates from various imputations
# active posture
mod_posture_mi <- with(data = imputed_data, exp = glm(active_posture ~ sense_power_factor + 
                                                        voted_putin + gender + age + education + income + rural_region,
                                                      family = "binomial"))
screenreg(combine_posture <- pool(mod_posture_mi))

# syria
mod_syria_mi <- with(data = imputed_data, exp = glm(syria_involvement_good ~ sense_power_factor + 
                                                      voted_putin + gender + age + education + income + rural_region,
                                                    family = "binomial")) 
screenreg(combine_syria <- pool(mod_syria_mi))

# crimea
mod_crimea_mi <- with(data = imputed_data, exp = glm(annexation_crimea_good ~ sense_power_factor + 
                                                       voted_putin + gender + age + education + income + rural_region,
                                                     family = "binomial")) 
screenreg(combine_crimea <- pool(mod_crimea_mi))

# --- ukraine
mod_ukraine_mi <- with(data = imputed_data, exp = glm(luhansk_donetsk_part_russia ~ sense_power_factor + 
                                                        voted_putin + gender + age + education + income + rural_region,
                                                      family = "binomial")) 
screenreg(combine_ukraine <- pool(mod_ukraine_mi))

# --- Table A14
screenreg(list(combine_posture, combine_syria, combine_crimea, combine_ukraine),
       booktabs = T,
       caption.above = T,
       caption = "OLS Results: Effect of Sense of Russian Power on Alternative Hawkish Measures",
       custom.coef.names = c("(Intercept)",
                             "Sense of Russian Power",
                             "Voted for Putin",
                             "Male",
                             "Age",
                             "Less than College",
                             "Income",
                             "Rural Resident"),
       stars = c(.001, .01, .05, .10),
       symbol = "\\wedge",
       custom.model.names = c("Assertive Posture in World Affairs",
                              "Military Involvement in Syria Effective",
                              "Annexation of Crimea Good for Russia",
                              "Luhansk and Donetsk Part of Russia within Ten Years"))


